function [R_maneuver,V_maneuver,t_maneuver,delV_total] = Run_StationKeeping_Linear(S_chase_eci,S_target_eci,S_target_ecef,...
    t_crossing,MAX_STATIONKEEPING_DURATION_HOUR,UTC_mission_start);
eps = 0.67;
mu = 398600.4418;
We = (2*pi / 86164.1) * [0;0;1];
options = odeset('RelTol' , 1E-6);
currentMinMin = 10000;
%% Get states 10 min apart for finding best delV in this duration
Eph = Ephemeris(S_chase_eci*1000, 200, 600,t_crossing);
Eph = Eph';
t = Eph(1,:);
R = Eph(2:4,:)/1000;
V = Eph(5:7,:)/1000;
%% Find the states of the chase vehicle in HCW frame of the S_target_eci;
% Define new variables
year = UTC_mission_start(1);
mon  = UTC_mission_start(2);
day  = UTC_mission_start(3);
hour = UTC_mission_start(4);
minute  = UTC_mission_start(5);
sec  = UTC_mission_start(6);
for ii = 1:length(Eph)
    UTC_new = [year, mon, day, hour, minute, sec+t(ii)];
    Rt=transpose(Eci2Ecef(UTC_new))*S_target_ecef(1:3);
    Vt=transpose(Eci2Ecef(UTC_new))*S_target_ecef(4:6) + cross(We,Rt);
    Rc = R(:,ii);
    Vc = V(:,ii);
    
    % Define HCW frame
    
    Ht = cross(Rt,Vt);
    %[km^2/s]Target vehicle specific angular momentum.
    
    Wt = Ht / dot(Rt,Rt);
    %[rad/s]Target vehicle angular velocity.
    
    Rhat = Rt / norm(Rt);
    %[]Current radial direction.
    
    Nhat = cross(Rt,Vt) / norm(cross(Rt,Vt));
    %[]Current normal direction.
    
    That = cross(Nhat,Rhat);
    %[]Current tangential direction.
    
    Hcw2Eci = [Rhat, That, Nhat];
    %[]Matrix that transforms vectors from HCW coordinates to ECI coordinates.
    
    Rct_rel = transpose(Hcw2Eci) * (Rc - Rt);
    %[km]Current relative position in HCW coordinates.
    
    Wt = transpose(Hcw2Eci) * Wt;
    %[rad/s]Transforms the target vehicle's rotational velocity from ECI to HCW coordinates.
    
    Vct_rel = transpose(Hcw2Eci) * (Vc - Vt - cross(Wt,Rct_rel));
    %[km/s]Current relative velocity in HCW coordinates.
    
    %% make initial point for finding smaller delV maneuver
    
    %% Determine best maneuver between given time interval
    
    MAX_STATIONKEEPING_DURATION_SEC = MAX_STATIONKEEPING_DURATION_HOUR*3600;
    
    n = sqrt(mu/norm(Rt)^3);
    
    deltat = 3600 : 100 : MAX_STATIONKEEPING_DURATION_SEC;
    totaldeltaV = zeros(1,length(deltat));
    deltaV1 = zeros(3,length(deltat));
    deltaV2 = zeros(3,length(deltat));
    for k=1:length(deltat)
        
        F = [4-3*cos(n*deltat(k)) , 0 , 0;...
            6*(sin(n*deltat(k)) - n*deltat(k)) , 1 , 0;...
            0 , 0 , cos(n*deltat(k))];
        %[] F matrix
        
        G = [(1/n)*sin(n*deltat(k)) , (2/n)*(1-cos(n*deltat(k))) , 0;...
            (2/n)*(cos(n*deltat(k)) - 1) , (1/n)*(4*sin(n*deltat(k)) - 3*n*deltat(k)) , 0;...
            0 , 0 , (1/n)*sin(n*deltat(k))];
        %[] G matrix
        
        Vct_rel_firstburn = G \ (-F * Rct_rel);
        %[] Velocity after the first burn
        
        deltaV1(:,k) = Vct_rel_firstburn - Vct_rel;
        %[] change in velocity 1
        
        Fdot = [3*n*sin(n*deltat(k)) , 0 , 0;...
            6*(n*cos(n*deltat(k))-n) , 0 , 0;...
            0 , 0 , -n*sin(n*deltat(k))];
        %[] F dot matrix. Derivative respect to deltat(k)
        
        Gdot = [cos(n*deltat(k)) , 2*sin(n*deltat(k)) , 0;...
            -2*sin(n*deltat(k)) , 4*cos(n*deltat(k))-3 , 0 ;...
            0 , 0 , cos(n*deltat(k))];
        %[] G dot matrix. Derivative respect to deltat(k)
        
        Vct_rel_secondburn = Fdot*Rct_rel + Gdot*(Vct_rel_firstburn);
        %[] Second burn at the meeting
        
        deltaV2(:,k) = -Vct_rel_secondburn;
        %[] delta V required for the second burn
        
        totaldeltaV(k) = abs(norm(deltaV1(:,k))) + abs(norm(deltaV2(:,k)));
    end
    
    delV = min(totaldeltaV);
    if delV < 0.003
        break
    end
end
min_index = find(totaldeltaV==delV);

Vct_rel_plus = Vct_rel+deltaV1(:,min_index);

deltat_min = deltat(min_index);

So_rel = [Rct_rel;Vct_rel_plus];

to = [t(ii):10:t(ii)+deltat_min];

[t_man , Sct_rel] = ode45(@(t,S)HCWModel(t,S,n) , to , So_rel,options);
t_man = t_man';
Sct_rel = Sct_rel';
Rct_rel = Sct_rel(1:3,:);
Vct_rel = Sct_rel(4:6,:);

[~,S_target] = ode45(@(t,S)R2bpCartesianModel(t,S,mu),to,[Rt,Vt],options);

%% Get all the states in Sct_rel from HCW to ECI frame.

S_target = S_target';
Rt = S_target(1:3,:);
Vt = S_target(4:6,:);
R_eci =zeros(3,length(Rt));
V_eci = zeros(3,length(Rt));
for k = 1:length(t_man)
    
    Ht = cross(Rt(:,k),Vt(:,k));
    %[km^2/s]Target vehicle specific angular momentum.
    
    Wt = Ht / dot(Rt(:,k),Rt(:,k));
    %[rad/s]Target vehicle angular velocity.
    
    Rhat = Rt(:,k) / norm(Rt(:,k));
    %[]Current radial direction.
    
    Nhat = cross(Rt(:,k),Vt(:,k)) / norm(cross(Rt(:,k),Vt(:,k)));
    %[]Current normal direction.
    
    That = cross(Nhat,Rhat);
    %[]Current tangential direction.
    
    Hcw2Eci = [Rhat, That, Nhat];
    %[]Matrix that transforms vectors from HCW coordinates to ECI coordinates.
    
    R_eci(:,k) = Hcw2Eci*Rct_rel(:,k)+Rt(:,k);
    
    Wt = transpose(Hcw2Eci) * Wt;
    %[rad/s]Transforms the target vehicle's rotational velocity from ECI to HCW coordinates.
    
    V_eci(:,k) = Hcw2Eci * Vct_rel(:,k) + Vt(:,k)+cross(Wt,Rct_rel(:,k));
    
end
n = length(t_man);
R_maneuver = [R(:,1:ii),R_eci(:,1:n-1),Rt(:,end)];
% V_maneuver = [V_eci(:,1:n-1),Vt(:,end)+0.01*rand];
V_maneuver = [V(:,1:ii),V_eci(:,1:n-1),Vt(:,end)];

t_maneuver = [t(1:ii),t_man];
delV_total = delV*eps;


end